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ABSTRACT 

It is widely believed that globular clusters evolve over many two-body relaxation times to- 
ward a state of energy equipartition, so that velocity dispersion scales with stellar mass as 
(T oc with rj = 0.5. We show here that this is incorrect, using a suite of direct N-body 
simulations with a variety of realistic IMFs and initial conditions. No simulated system ever 
reaches a state close to equipartition. Near the center, the luminous main-sequence stars reach 
a maximum r/max ~ 0.15 ± 0.03. At large times, all radial bins convergence on an asymptotic 
value ?7oo ~ 0.08 ± 0.02. The development of this "partial equipartition" is strikingly similar 
across our simulations, despite the range of different initial conditions employed. Compact 
remnants tend to have higher 77 than main-sequence stars (but still 77 < 0.5), due to their 
steeper (evolved) mass function. The presence of an intermediate-mass black hole (IMBH) 
decreases rj, consistent with our previous findings of a quenching of mass segregation under 
these conditions. All these results can be understood as a consequence of the Spitzer instabil- 
ity for two-component systems, extended by Vishniac to a continuous mass spectrum. Mass 
segregation (the tendency of heavier stars to sink toward the core) has often been studied 
observationally, but energy equipartition has not. Due to the advent of high-quality proper 
motion datasets from the Hubble Space Telescope, it is now possible to measure 77 for real 
clusters. Detailed data-model comparisons open up a new observational window on globular 
cluster dynamics, structure, evolution, initial conditions, and possible IMBHs. A first compar- 
ison of our simulations to observations of Omega Cen yields good agreement, supporting the 
view that globular clusters are not generally in energy equipartition. Modeling techniques that 
assume equipartition by construction (e.g., multi-mass Michie-King models) are approximate 
at best. 

Key words: stellar dynamics — globular clusters: general — methods: n-body simulations 



1 INTRODUCTION 

Gravitational encounters within stellar systems in dynamical equi- 
librium, such as globular clusters, drive evolution over the two- 
body relaxation timescale. The evolution is toward a thermal ve- 
locity d istribution, in which stars of different mass have the same 
energy (Spitzeij l 19871) . This thermalization also induces mass seg- 
regation. As the system evolves toward energy equipartition, high 
mass stars lose energy, decrease their velocity dispersion and tend 
to sink toward the central regions. The opposite happens for low 
mass stars, which gain kinetic energy, tend to migrate toward the 
outer parts of the system, and preferentially escape the system in 

the presence of a tidal field. 

As first pointed out by ISpitzetl d 1969I) . not all self-gravitating 
systems can attain complete energy equipartition, i.e., a mass- 
dependent velocity dispersion scaling as a{m) oc m"'^'^. For ex- 
ample, in a simple two-component model with light and heavy stars 
of masses mi and m2, respectively, equipartition is possible only if 



E-mail addresses: trentiast.cam.ac.uk 



the mass fraction in heavy particles is smaller than a critical value; 
M2/M1 < 0.16(mi/m2)^/^, where Ah and M2 are the total 
masses in stars of mass mi and m2, respectively (IS pitzerl [19691) . 
With too many massive particles, no equipartition is possible, and 
the heavy particles tend to become a self-gravitating system, pro- 
ducing a deep core collapse while continuing to transfer energy to 
the lighter stars. 

Subsequent research has refined our knowledge of the con- 
ditions required to establish energy equipartition. This was stud- 
ied for two-component systems using a variety of numerical 
met hods, including pio n eering work by ISpitzer&Har3 ( Il97ll) 
and Inaga ki & Wivantg jl984h with Fokker-Planck and Monte 
Carlo schemes, followe d by more realistic simulations both with 
a Monte Carlo c ode jWatters e t al. 200 0}) and with direct N- 
body integration jKhaUsi et 31.112007.) . .Vishniad ( Il978h general- 
ized the analytical Spitzer instability analysis to include a contin- 
uous mass spectrum. He showed that cl usters with realistic mass 
functions, such those of the Salpete] ( Il955h form, are unable 
to attain energy equipartition. With a complementary approach 
based on multi-component King distribution function modeling, 
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iKondrat'ev & Ozemovl [m^ showed that even if a system attains 
equipartition in the energy phase space, the mass-dependent kine- 
matics does not necessarily scale as a{m) oc m~^'^. 

Much of the theoretical interest in equipartition has focused 
on how energy exchange affects the process and time scale of core 
collapse in globular clusters. Both analytically and numerically, 
it has been established that the core collapse timescale in a two- 
component model is inversely proportional to the ratio of the heavy 
to light particle m asses (e.g.. IPortegies-Zwart & McMillanlbOOd : 
iKhalisi et alj2007l ). This has the important consequence that young 
star clusters can undergo a runaway central collapse of massive 
stars within the first few million years after their formation, pos- 
sibly providing a pat hway to the formation of i ntermediate-mass 
black holes (IMBHs) jPortgies-Zwart et alJZOoi) . 

One of the directly observable consequences of the drive to- 
wards energy equipartition is spatial mass segregation. In simula- 
tions, the amount of mass segregation reaches a steady-state con - 
figuration within a few initial relaxation times ( lGiUetal.ll20"08h . 
The presence of an IMBH, or of a stellar-mass black hole binary, in 
the core of a s tellar cluster reduces the amount of mass segregation 
that develops ((Saumgardt et alj|2004l : iTrenti et al ] |2007tlGillerai] 
l2008l : IUmbreit & RasidboiiT 

Mass segregation has been measured in several globular clus- 
ters using high-quality Hubble Space Telescope (HST) data. This 
only requires observations of stellar positions and luminosities (on 
the main sequence, luminosity correlates with mass). Such obser- 
vations have confirmed the qualitativ e picture that massive stars 
are preferentially found in the core (|de Marchi & Parescel Il994l : 
lAndreuzzi et al.ll200(]|: [Paust et al.ll201oir However, interpretation 
of such data requires detailed dynamical models. These are of- 
ten constructed assuming energy e quipartition for all components 
throiigh Michie- King mode ls (e.g.. lGunn & G riffin 1979; Mevlan 
Il988| - [de Marchi et al.ll2000 '). The accuracy of this is unclear, given 
the theoretical results described above which indicate that globular 
clusters may not generally be in energy equipartition. 

A preferred approach is to compare data directly to the re- 
sults of numerical simulations. Detailed, data-model comparisons 
for globular clusters such as NGC2298 and MIO have shown a re- 
markable quantitative agreeme nt with the mass segregation predic- 
tions f rom N-body simulations jPasguato et al.l200^lBeccari et al.l 
l20ld : lUmbreit & Rasid I20TI . However, such comparisons are 
challenging to carry out for general datasets. It not only requires de- 
tailed numerical simulations, but the simulations must also be "ob- 
served" to mimic the magnitude limits, completeness and crowding 
effects in the data. These observational effects generally vary with 
radial distance within a system, and between systems. 

Mass segregation provides only an indirect measure of en- 
ergy equipartition: stellar energies can only be estimated statisti- 
cally from stellar positions, and this requires assumed knowledge 
about the dynamical state and gravitational potential of the system. 
It has not been possible historically to measure a{m) in globu- 
lar clusters directly from knowledge of individual stellar velocities. 
Line-of-sight velocities of stars can be measured from spectra, but 
those are difficult to obtain for faint stars below the main-sequence 
turn-off. Stars in globular clusters with measured line-of-sight ve- 
locities are therefore generally giant stars. These stars all have very 
similar mass (because stellar evolution proceeds rapidly after the 
main-sequence turn-off), making it impossible to determine how 
the velocity dispersion cr in a cluster depends on stellar mass m. 

This situation is now changing due to the advent of HST stud- 
ies of stellar proper motions in globular clusters, which can quan- 
tify the kinematics of the stars as function of their position along the 



main sequence (i.e., their mass). [Anderson & van derMarellfcOld) 
presented measurements for more than 100,000 stars in the central 
region of Omega Gen, finding cr oc m^'^'^ . Proper motion cata- 
logs for some two d ozen additional clusters are under constructions 
( lBellinietalJl2013l) . This opens up a new observational window 
into the dynamics of globular clusters. 

To understand and exploit these new data there is a need for 
new models that predict in detail the mass- and spatially-dependent 
kinematics of globular clusters as a function of initial conditions 
and evolutionary state. With new data and models in hand, there 
is the potential to validate our understanding of globular cluster 
structure and evolution at a new level. Moreover, knowledge of the 
cr(m) relation may constrain structural parameters that are diffi- 
cult to assess otherwise, such as the possible presence of a central 
IMBH. 

The goal of this paper is to make the first step towards mod- 
eling mass-dependent kinematics in the era of detailed proper mo- 
tions measurements. We analyze a sample of realistic direct N-body 
simulations of star cluster evolution, and we also use the models 
to interpret the specific observational results obtained for Omega 
Cen. The structure of the paper is as follows. In Section |2] we in- 
troduce the numerical simulations that we have carried out, and we 
describe how we analyze simulation snapshots to construct cr(m). 
The results of our analysis are presented in Section |3] We compare 
the model predictions to Omega Cen observations in Section^ and 
find excellent consistency between models and observations. We 
summarize our key findings and conclude with an outlook for fu- 
ture studies in Section|5] 



2 N-BODY MODELS 

2.1 Setup and Methodology 

We analyze the stellar d ynamics in the dire ct N-body simulations 
previously carried out bv lTrentietai1(l20I0l) to investigate the evo- 
lution of structural parameters in star clusters. For full details on 
the code and the simulation setup we refer to that paper (and previ- 
ous investigations in t he s ame framework by [ Trenti. Heggie & Hull 
l200llGill et al.l2008l. and lPasguato et al.l2009l) . To summarize, we 
follow the dynamical evolution of simulated star clusters using the 
NB0DY6 code tAarseth.2003t) , which guarantees an exact treat- 
ment of multiple interactions between stars by employing special 
regularization techniques, without resorting to the introduction of 
softening. The dynamics of the system are thus followed with ex- 
tremely high accuracy, at the price of a using lower number of par- 
ticles compared to more approximate methods (e.g., Monte Carlo 
codes). 

We simulate systems with TV = 32768 to A'^ = 65536 parti- 
cles, with an initial mass function constructed to reproduce the long 
te rm evolution of star clusters that are > 10 Gyr old. As desc ribed 
in iTrenti et al.l ( l2010t) we start from either a lSalpeteij(ll955h or a 
[Miller & Sc^ ( Il979l) initial mass function a nd then apply an in- 
stantaneous step of stellar evolution with the [Hurley et al.l i2000l) 
evolutionary tracks to obtain a main-sequence tumoff at 0.8 M©. 
Our standard assumption is to retain 100% of the dark remnants 
without assigning them velocity kicks, but we explore a lower re- 
tention fraction (30%) of neutron stars and stellar mass black holes 
in two runs. During the simulation, only gravitational interactions 
are considered and stars are not evolved further. This simplified ap- 
proach is ideal to highlight the effects of gravitational forces on 
the long-term evolution of star clusters and the drive toward energy 
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Table 1. Summary of N-body simulations and energy equipartition results 
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N-body simulations of star clusters in a tidal field with self consistent King model initial condi- 
tions. Column (1): number of particles A'. Column (2): binary fraction /. Column (3): initial mass 
function used (Sa: Salpeter power law; MS: Miller & Scalo). Column (4): initial concentration 
of the density profile (King index Wq)', Columns (5)-(8): equipartition results for single main- 
sequence stars, quantified as desciibed in Section l3.1.2l Column (5): time tjiaif to reach half of the 
maximum equipartition; Column (6): time tmax to reach the maximum equipartition; Column (7): 
maximum equipartition r^max reached, where a{m) oc . Column (8): late-time asymptotic 
equipartition r)ao of the system as a whole. Columns (5)-(7) pertain to the innermost 10% of the 
stars as seen in projection, while column (8) pertains to the system as a whole. 

Table notes. *: canonical simulation discussed in Section [3l1 f: 30% NS/BH retention; |: 
"^IMBH = O.Olrriciustor; §: "/"t = 6.28 [compact initial conditions where the Roche lobe is 
under-filled by a factor 2]. 



equipartition, without the added complication of disentangling dy- 
namics from energy injection induced by stellar evolution (through 
loss of mass with negative energy). Our approach provides an ap- 
proximate, yet accurate description of the current evolution of glob- 
ular clusters, since stellar evolution affects star cluster structu re sig - 
nificantly only during the first few hundred Myr (e.g. see iHurlevi 
l2007l;lMackev et alJlOOSi) . 

The initial conditions in the position and velocity space are 
drawn from a King distribution function (King 1966), with scaled 
central potential Wo ~ 3, 5, 7. Initially, all particles of mass m 
are drawn from the same velocity distribution, hence at any given 
radius a{m) is constant. Our runs include a primordial binary frac- 
tion / — Nt/{Ns + Nt) ranging between and 0.1, with Ns and 
Nb being the number of singles and binaries respectively. All bina- 
ries are "hard" in the lHeggi3 ( Il975t) classification (i.e., semi-axis 
typically smaller than ~ 10 AU for a typical globular cluster with 
central velocity dispersion of 10 km/s). One of our A*' — 32768 
runs contains a central IMBH, with the BH mass set at 1 % of the 
total cluster mass. Table[T]gives an overview of all the simulations 
performed and analyzed. 

The simulated star clusters are tidally limited and particles ex- 
perience a tidal force from a point-like parent galaxy, assuming 
that the cluster is in circular orbit at a distance selected so that the 
tidal radius is self-consistently defined by the King density profile. 
Hence, models have different tidal field strengths for different Wo 



values (see lTrenti. Heggie & Hij? '2007 and lTrenti et"aLll2010l for a 
full definition of the tidal field equations). 

The simulations are run for t > 15 trh(O), where trh{0) is 
the initial half-mass relaxation time, but they are stopped earlier 
when 80% or more of the initial mass in the system has been lost 
due to evaporati on of stars. In the natural (dimensionless) N-body 
units defined bv lHeggie & Mathieij J 19861) the relaxation time can 
by written as: 

_ai38iVr^ 

" logiO.llN) ' ^ ' 

where rn ~ 1 is the half-mass radius and the time unit approx- 
imately corresponds to the orbital period of a particle at r^. For 
simplicity, all our time-dependent results are expressed in units of 
this relaxation time, since this is the relevant timescale for develop- 
ment of energy equipartition and thermodynamic equilibrium. 

2.2 Structural Evolution 

We refer to lTrenti et al.l(l2010l) for a detailed discussion of the evo- 
lution of the cluster structural parameters. To summarize, after a 
few relaxation times, all clusters tend to evolve toward a quasi- 
universal core-halo structure with a central concentration (mea- 
sured either as the core to half mass radius ratio rc/r^ or as core 
to tidal radius ratio Vc/rt) which is independent of its value at the 
beginning of the run. The long-term equilibrium value for the con- 
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centration is defined uniquely by the efficiency of the central energy 
production due to dynamical processes: runs with only single parti- 
cles undergo a marked core collapse and have higher concentration, 
while in tuns with primordial binaries the core collapse is halted at 
lower central densities because of efficient energy input by three 
body encounters. The run with an IMBH has the largest core at late 
times. Once the long term Vc/rh ratio has been established, the sys- 
tem undergoes a self-similar expansion at fixed concentration, until 
the majority of the stars have been lost by tidal evaporation, or the 
mechanism of central energy production is depleted (for example 
because of disruption of all primordial binaries). 

This classical picture for the dynamical evolution of a 
globular cluster is clearly present when one considers three- 
dime nsional, mass-based definitions of the core and hal f-mass ra- 
dius jTrenti, Heggie & Hut 2007l: lFregeau & Rasiol2007h . The situ- 
ation becomes more complicated when simulations are instead "ob- 
served", namely when only two-dimensional, light-based defini - 
tions for structural parameters are constructed (tTrenti et alj|2010l) . 
For example, it is then possible that clusters without primordial bi- 
naries still show a large core at late times if the core is dominated by 
heavier dark remnants (neutron stars and stellar mass black holes). 



2.3 Analysis of Stellar Dynamics 

Our analysis here differs from that in lTrenti et a l. in that we 

now consider the stellar dynamics in the simulations. Simulation 
snapshots with particle positions and velocities were saved every 
10 dimensionless time units, yielding several tens of snapshots per 
relaxation time. For each snapshot, at time t, we binned the par- 
ticles to derive the velocity dispersion a{m,r,t) as a function of 
mass and projected radius. Mass binning was carried out in inter- 
vals of 0.1 M0. Two-dimensional (projected) radial bins were de- 
fined by Lagrangian radii, each containing 10% of the particles in 
the snapshot. 

We took the following steps to minimize noise from low par- 
ticle numbers. First, we measured three-dimensional velocity dis- 
persions, and then rescaled the results to two dimensions with a 
multiplying factor 1^/2/30 Second, we created projections for the 
radial distribution of the stars along three orthogonal directions and 
averaged them. And third, we averaged together 30 independent 
neighboring time snapshots. These steps allow us to achieve an ef- 
fective number of particle velocities in our analyses that for our 
A'^ = 65536 runs is comparable to the number of stars in a massive 
globular cluster like Omega Cen. 

We determined the functional dependences of a(m, r, t) sep- 
arately for single main-sequence stars, binary stars, and compact 
remnants (white dwarfs, neutron stars and stellar mass black holes) 
in the simulations. For comparison to observational results, we 
focus in the following primarily on the results for single main- 
sequence stars. In reality, some stars in observational catalogs may 
be unresolved binaries. However, we found that our main results 
described below did not change significantly if binary stars in the 
simulations were included in the single-star analyses (i.e., treating 
each binary of total luminosity L as a single star of that luminosity). 



This procedure implicitly assumes isotropy in the three-dimensional ve- 
locity distribution, as expected for a collisionally relaxed stellar system. We 
did verify that our results do not change (except for somewhat increased 
noise) if we compute two-dimensional velocity dispersions directly. 



3 EQUIPARTITION RESULTS 
3.1 Canonical Simulation 

We describe first the results for a "canonical" simulation: a model 
with = 65536 stars, a Miller & Scalo initial mass function 
(IMF), no primordial binaries, and a King concentration parameter 
Wo = 5 (see Table [TJ. This model has a moderate ini tial concen- 
tration and external tidal field strength. As discussed in lTrenti et al.l 
( 2010), the structural evolution of this cluster is fairly proto-typical. 
Starting from an initial concentration c = — logjQ(rt/rc) ~ 1, it 
undergoes core-collapse in about t ~ 7trh(0). At this stage, its ob- 
served concentration (based on the core radius observed from the 
surface brightness profile) has reached a quasi-equilibrium value 
c ~ 1.7. The three-dimensional concentration is instead higher 
(c > 2.5) as the 3D core contracted to reach densities suffi- 
cient to produce a few binaries through three body encounters. 
At t > ItrhiQ), the system continues to evolve with self-similar 
density and surface-brightness profiles, progressively losing mass 
as a result of the external tidal field until complete dissolution at 
f ~ 16t,h(0). 



3.1.1 Mass Dependence of Kinematics 

Figure [T] shows a{m), determined as described in Section 12.31 
for the innermost radial Lagrangian bin (10% enclosed mass) at 
t = 5.1 trh{0). This is during the initial core contraction as the sys- 
tem is on its way toward core collapse. Red points pertain to single 
stars along the main sequence, while blue points pertain to compact 
remnants which are generally heavier. Single stars are found up to 
the tumoff mass of 0.8 Mq, while the lightest compact objects have 
masses of ~ 0.55 Mq. 

The main sequence stars and the compact remnants have com- 
parable dispersions over the mass range where both types of objects 
exist, 0.55 < m < 0.8 Mq. However, when the full range is con- 
sidered over which each object type is found, the compact remnants 
have a significantly steeper a{m) relation than the main-sequence 
stars (approximately m-o-^* vs. m-«-", respectively). Neither 
type of object achieves energy equipartition, but the heavier com- 
pact objects are closer to it than the lighter main-sequence stars. 
This behavior is qualitatively consistent with the stability analysis 
presented by Vishniac (1978) , given that the remnants have a sig- 
nificantly steeper mass function than the main-sequence stars (see 
Figure|2j. This makes it easier for the population in the 0.7 — 2 Mq 
mass range to thermalize. The steep mass function for remnants is 
due to the significant mass loss experienced by massive stars before 
they become remnants, and is relatively insensitive to the adopted 
IMF. The generic shape of a{m) in Figure[T]is therefore typical for 
all our simulations. 

The relation a{m) in Figure[T]shows a break at m ~ 0.7 Mq. 
Figure |2] shows that this is also the mass at which the density of 
main-sequence stars is equal to the density of the remnants. The 
IVishniacf ( Il978l) analysis suggests that more energy equipartition 
is possible for a sub-population within a given mass range, when 
the lighter counterparts are much more abundant. So thermaliza- 
tion is more likely for objects above 0.7 M0, as there are plenty of 
lower mass particles along a steep mass function {m~^'^). By con- 
trast, the lowest mass remnants and the main sequence stars below 
0.7 M0 are less able to thermalize efficiently, because the mass 
function is only increasing mildly (m~^'^). This explains the quan- 
titative results in Figure [T| 

Figure [T] shows that the a{m) relations for single main- 
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Figure 1. Velocity dispersion ct as a function of particle mass m around 
time t = 5.1trj,(0) for particles at the center of the system (two- 
dimensional projected radius r < 0.78 r^) in our canonical A'^-body sim- 
ulation (the N = 65536 nin with Wo = 5, a lMiller & Scalol {\9l4i IMF, 
and no primordial binaries). The red points with errorbars are for single 
main sequence stars, while the blue points are for compact remnants. Lines 
show the best fitting power-law relation cr{rn) oc m^^ for each set of ob- 
jects, with r} = 0.14 for single main sequence stars and r} = 0.34 for 
compact remnants. 
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Figure 2. Mass function for the particles of an A'^-body run (including, e.g., 
our canonical simulation) starting from a Miller & Scalo 1979 IMF, evolved 
to a turn-off mass of 0.8 Mq. The red histogram represents main-sequence 
stars, and the blue one compact remnants. Approximate power-law fits are 
show as dashed curves, which have slopes m~^'^ and m~^ '^, respectively. 
Since the mass function for remnants is much steeper, it is easier for them 
to approach energy equipartition (see Figure[T]. 



sequence stars and compact remnants separately, are both reason- 
ably well fit by a power law of the form (j{m) oc m"''. The 
compact remnants are only of theoretical interest as they are not 
directly observable. In the following, we therefore focus our atten- 
tion only on the main sequence stars, which can be observed. Using 
the a{m, r, t) profiles constructed as described in Section [23] we 
fit for each radial Lagrangian bin r and at each time t the mass- 
dependent kinematics assuming a power law cr(m) cx over 
the mass range 0.2 m/MQ ^ 0.7. Here is a free parameter, 
with Tj = Q corresponding to absence of equipartition (velocity dis- 
persion independent of mass), Q < rj < 0.5 corresponding to par- 
tial equipartition, and 77 = 0.5 corresponding to complete equipar- 
tition. The quantity is a well-defined fit-parameter that can be 
compared across models, even when the a{m) relation shows some 
residual curvature, as e.g. in Figure [T] 

3.1.2 Radial and Time Dependence of Equipartition 

The degree of energy equipartition as a function of time, rj{t), is 
shown for different Lagrangian radii in the canonical simulation in 
Figure |3] At t = 0, the system starts (by construction) with no 
equipartition at all. As time progresses, there is a relatively rapid 
linear rise in rjit) for the inner Lagrangian radii. This is followed 
by a flattening toward a maximum value around 77 ~ 0.15 that is 
reached after a few initial half-mass relaxation times. Subsequently, 
the value of rj drops slowljQ. For the outer Lagrangian radii, the 
value of rjit) evolves on a longer timescale, reflecting the longer 
relaxation time in the outskirts of the system. For all radii there is 
convergence to a value around 77 ~ 0.10 at late times. However, be- 
cause of the slower development of equipartition, the outer regions 
do not first "overshoot" this value by reaching a maximum at early 
times, as seen for the inner Lagrangian radii. 

To describe the time evolution of equipartition empirically, we 
adopted a fit function of the form: 



where 



Vtitit,r) = r]c{r) x fa[t/tcir)], 



(2) 



(3) 



Here tc{r) and ric{r) are a characteristic time and a characteristic 
rj value, both of which depend on radius. The parameter a{r) also 
depends on radius and defines the exact shape of the function fa. 
Upon varying tc, r]c and a for each radial bin to optimize the fit, 
this empirical formula generally provides an adequate description 
of the time dependence of equipartition in the simulations. For the 
canonical simulation, this is illustrated by the solid curves in Fig- 
ure [3] 

To compare the results of different simulations, it proved 
convenient to extract the following quantities, which capture the 
essence of the time evolution of equipartition. The quantities were 
determined for all simulations from the empirical fits of the form 
given by equation (2] rather than from the individual simulation 
snapshot datapoints, because the latter are available only at discrete 
times. 

• ??max: the maximum rj that is reached in the central region, 
defined here as the innermost Lagrangian radius (inner 10% of the 
projected mass). 

^ This drop in r] is likely due to formation of binaries in the core. These 
provide heating that acts to reduce mass segregation and the approach to 
equipartition. 
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Figure 3. Time evolution of the energy equipartition indicator rj (the power- 
law slope in the relation between velocity dispersion and mass, a oc m^^) 
for single main sequence stars in our canonical Af-body simulation. Data 
points with error bars were determined directly from the simulation snap- 
shots, as described in Sections |2.3| and 13.1.11 SoHd curves are fits of the 
form given by equation l|2). Red refers to the innermost 10% of the main 
sequence stars in projection; blue refers to three projected annuli contain- 
ing 10% of the particles each, and ending at 30%, 50% and 70% of the 
enclosed particle number; green refers to the outer Lagrangian bin at 90% 
of the enclosed particle number. The time along the abscissa is expressed 
in units of the initial half-mass relaxation time t^f^ (0). These results for the 
canonical simulation are representative of the general behavior of our N- 
body simulations. The inner region of the system reaches a higher degree 
of equipartition within a few relaxation times than the outer regions. An in- 
crease in Tj takes longer to develop at large radii. However, at later times, 
the whole system converges toward a common value of r). Arrows indicate 
the characteristic times defined in Section l3.1.2l Complete energy equipar- 
tition (r) = 0.5) is never attained, confirming previous investigations based 
on stability analysis. 

• tmax: the time at which the central region achieves r] = r^max. 

• thau: the time at which the central region achieves 77 = 
0.5r;max. This value is defined more robustly than tmax, since rj{t) 
tends to flatten near its maximum. 

• too '■ the late time t > tmax at which the 77 values for all La- 
grangian radii are most similar (defined as minimum dispersion in 

n)- 

• rjao'- the value of 77 at t = too averaged over all Lagrangian 
bins. This represents approximately the long-term asymptotic value 
of T) for the system as a whole. 

We list these quantities in Table [T] for all simulations. We omit 
too since it is not very robustly determined, lacks a clear physi- 
cal meaning, and may depend on the (arbitrary) time at which the 
simulation is ended. Since ri{t) is roughly flat at large times, 7700 by 
contrast is defined more robustly. 

For the canonical simulation, thaif = 0.9trh(0), tmax = 
<iAtrh{Q), and too = 12.5f,;i(0). These times are indicated with 
arrows in Figure |3] Half of the maximum equipartition develops 
rather quickly in the central region of the system, while the maxi- 



Figure 4. Scatter plot of equipartition results, showing Tjmax vs. thalf for 
all of the Af-body simulations. The velocity dispersion of single main- 
sequence stars scales as cr oc m^^, and the quantity r/max is the maxi- 
mum power-law slope r] for the innermost 10% radial Lagrangian bin in 
projection. The quantity t^aif is the time at which this radial bin reaches 
77 = 0.57)max. Simulations with N = 32k particles are shown in red, and 
those with = 64fc particles in blue. The initial rise in 77(4) is approxi- 
mately Hnear, so all simulations have essentially achieved a near-maximum 
equipartition in their cores after t > 3 trh{0)- 

mum value is reached only around half-way to core-collapse. The 
maximum equipartition is characterized by 77niax ~ 0.144 ±0.006, 
and the long-term asymptotic value is 7700 — 0.085 ± 0.006. 

The analysis discussed above was carried out using pro- 
jected radii, which are the only observable quantities. The results 
show that the canonical simulation never attains complete energy 
equipartition. However, at early times the particles in the projected 
core are closer to complete equipartition than what is found for 
all particles at late times. So it expected that stars in the three- 
dimensional core might achieve an even higher degree of equipar- 
tition at early times. We analyzed rj also for the three dimensional 
core, and found that this is indeed the case. However the effect 
is very modest, and the three-dimensional core is only marginally 
closer to equipartition (by A77 ~ 0.05 at t = tmax)- Therefore, no 
region in the simulation is attaining complete energy equi partition, 
in full agreement with the analytical stability analysis of Ivishnia3 
( fr978l) . 

3.2 Dependence on Model Parameters 

3.2.1 Sample Statistics 

We have found that the analysis of energy equipartition across our 
sample of A'^-body simulations shows a surprisingly uniform and 
consistent picture. All the runs behave to first approximation like 
the canonical run discussed in Section [3T| Figures |4] and |5] show 
scatter-plot representations the results for thaif, tmax, 77max and 
7700 for all of the simulations. The uncertainties in rimax and r^oo 
are generally small, A77 < 0.01, and are smaller than the scatter 
between different simulations. 
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Figure 5. Scatter plot of equipartition results, showing rjao vs. imax for 
all of the A'^-body simulations. The velocity dispersion of single main- 
sequence stars scales as it oc m~^, and the quantity rjoo is the asymp- 
totic power-law slope toward which the system as a whole converges at 
late times. The quantity tmax is the time at which the innermost 10% 
radial Lagrangian bin in projection achieves its maximum value rjmax. 
Simulations with A'^ = 32k particles are shown in red, and those with 
N = 64fc particles in blue. The quantity r?oo has a near-universal value of 
0.08 ± 0.02 across all simulations, and is very far from complete equipar- 
tition (»7 = 0.5). 



The value of ?7max varies from rjmax = 0.12 achieved by 
the run with a central IMBH and by a run with low concentration 
(Wo = 3) and primordial binaries (/ = 0.055) to rjmax = 0.18 
for the runs with a low retention fraction of neutron stars and stel- 
lar mass black holes. Runs with a higher fraction of primordial 
binaries develop less equipartition than their similar counterparts 
with single stars only, although the effect is overall modest (e.g. 
from T]max = 0.159 to rjmax = 0.127 going from / = to 
/ = 0.1 for the serie s of ru ns with iV = 32768 particles, Wp = 7 
and a Miller & S cal9 (^1979) IMF). Comparing runs with a lSalpetej 
versus Ivliller & Scalo ( 1979) IMF does not show any sig- 
nificant difference. Table[T]also shows that lower initial concentra- 
tion Wo generally yields slightly lower r]max ■ As discussed in Sec- 
tion [3]2]2] below, this is most likely due to the different tidal field 
strength assumed in simulations with different Wo, rather than due 
to the concentration itself. 

The long-term asymptotic value r^oc is even more uniform 
across the sample of simulations than is 77max- Its value ranges from 
rjoo = 0.063 to rjoo ~ 0.101. The sensitivity of 7700 to model pa- 
rameters is qualitatively similar to that for ?7max, and there is a mild 
correlation between the two quantities across the sample. A Spear- 
man rank correlation test gives a correlation r = 0.36, which is 
significant at greater than 90% confidence for our sample of 18 
simulations. Since some of the mechanisms that drive the suppres- 
sion of the equipartition are most active in the core rather than in the 
whole system, such as the presence of an IMBH, it is not surprising 
to observe only a modest correlation. 

The timescale needed to approach maximum equipartition, as 



measured by either fhaif or fmax, varies more between simula- 
tions. The value of fhaif/irh(0) ranges from approximately 0.5- 
1.6. The initial rise in rj{t) is approximately linear, so this im- 
plies that after t > 3 trh{0) all systems have essentially achieved 
a near-maximum equipartition in their cores. The actual values 
of fmax/trh(0) show Considerable variation though, ranging from 
2.1-7.2. This is because the ri{t) profile beyond the maximum is 
often nearly flat for the innermost Lagrangian bin, so that tmax is 
not robustly determined. An example of this is provided in Fig- 
ure |6l which is similar to our canonical simulation in Section [3T1 
but pertains to an initial King concentration index Wo ~ 7. 

Figure |4] shows that the runs with A'^ — 65536 particles tend 
to have longer fhaif/tr/i(0) than those with = 32768. However, 
this may be due to small systematic errors in the determination of 
ihaif- Figure |6] shows that the steep initial rise in rj{t) in the inner 
Lagrangian bin is not always well reproduced by the fit of the func- 
tional form in equation (|2}- The values of r/max, rjoo and tmax do not 
show any clear trends with particle number. This suggests that our 
simulations have sufficient A'^ to determine real physical trends (as 
opposed to numerical artifacts, transients, or shot-noise dominated 
features), and that any dependence on A^ is properly accounted for 
by expressing time scales in units of trh{0). 

The equipartition timescales fhaif and tmax do not show ob- 
vious correlations with other parameters or with the model initial 
conditions. There is a hint of negative correlation between fhaif and 
T?max (Figure|4}, but we are hesitant to attach much significance to 
this, given the possibility of small systematic errors in the determi- 
nation of thaif • No obvious correlation is present between tmax and 
rioo (Figure|5}. 

3.2.2 Interpretation 

The maximum equipartition value T^max correlates with model ini- 
tial conditions in a way that can be understood based on the phys- 
ical insights, stability analyses, and mass segregation research dis- 
cussed in Section[T] 

It has been shown previously that the presence of an IMBH 
suppresses mass segregation. So the effect is naturally to reduce 
the amount of energy equipartition, measured by r), as well. Fig- 
ure|7]shows the rjit, r) profiles for our simulation with an IMBH. 
These can be compared to the results in Figure|6]for a similar model 
(with twice the number of particles) without an IMBH. The inclu- 
sion of an IMBH yields a modest reduction in ?7max- We attribute 
this to the fact that the IMBH generally has at least one particle 
tightly bound to it (often a compact remnant). This binary system 
scatters single main sequence stars out of the core, independently 
of their mass (as these third bodies have masses much smaller than 
the IMBH binary). This partially counters the natural tendency for 
more massive main-sequence stars to segregate into the core and 
have a lower velocity dispersion compared to lighter counterparts. 

Primordial binaries also act toward suppressing mass segrega- 
tion, with a qualitatively similar mechanism. So it is not suipris- 
ingly that we generally find a decrease in rimax with increasing / 
(see Table [T](. And stellar mass black holes tend to form BH-BH 
binaries, which also act in a similar way. Hence, it is not surprising 
that our two simulations with a low retention fraction of neutron 
stars and BHs, yield the highest values of rymax (i.e., achieving the 
most equipaitition in their central region). 

Finally, we find that runs with a steeper IMF yield similar 
T/max b ut marginally hi gher 7700- This is consistent with the anal- 
ysis of lVishniad h978l) . With a steeper IMF, heavy particles can 
transfer their kinetic energy more efficiently to their lighter coun- 
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Figure 6. Time evolution of the energy equipartition indicator r\ (the power- 
law slope in the relation between velocity dispersion and mass, cr cx: m^'') 
for single m ain sequence sta rs in a simulation with N = 65536 particles, 
Wq = 7, a lSalpete3 J 19551) IMF, and no primordial binaries. The overall 
evolution is similar to that for the canonical simulation (see Figure|3] which 
uses the same layout and color-coding). However, the present siinulation 
shows "energy divergence" in the outer parts at early times, with evolution 
toward r] < (i.e., high-mass stars have higher velocity dispersion than 
low-mass stars). This transient behavior is possibly related to ejection of 
particles from the core into the outer parts of the system because of the 
dynamical interactions (see Section |331 . This simulation also shows that 
the best-fit function of the form given by equation l|2) may underestimate a 
very rapid rise in r](t) at early times, yielding a value of t^aif that is biased 
high. 

terparts, since those are more numerous. Hence, more equipartition 
can be reached. 

The theoretical expectation is that the initial concentration 
of the system should not play a significant role in determining 
This is because after several relaxation times the system 
has evolved toward a quasi-equilibrium state with a u niversal con- 
centra tion, independently from where it started from dTrenti et al.l 
I2OIOI) . As mentioned in Section [3.2.1| our results taken at face value 
do not confirm this. We find lower r^max in simulations started with 
lower Wo- However, this is most likely because of the effects of 
the tidal field on the kinematics of the system. The concentration in 
the initial conditions of most of our simulations is self-consistently 
associated with the assumed tidal field strength. Instead, our simu- 
lation starting from low initial concentration Wo = 3, but with an 
under-filled tidal Roche lobe, yields similar r^max as a simulation 
starting with Wo = 7. So we conclude tentatively that concentra- 
tion does not affect equipartition, while a strong tidal field has a 
mild effect of suppressing the maximum rj. 

3.3 Energy Divergence 

The expectation based on generic thermodynamic arguments is that 
clusters evolve toward energy equipartition (even if they may never 
reach complete equipartition). So if one stars from a state with 
7] = (velocity dispersion independent of mass m), one expects 



Figure 7. Time evolution of the energy equipartition indicator rj (the power- 
law slope in the relation between velocity dispersion and mass, cr oc m~'') 
for single main sequence st ars in a simulation wit h a central IMBH, N = 
32769 particles. Wo = 7, a lMiller & Scalol i 1979 1) IMF, and no primordial 
binaries. The overall evolution is similar to that for the corresponding sim- 
ulation without a central IMBH (see Figure[6] which uses the same layout 
and color-coding). However, the IMBH suppresses the ainount of equipar- 
tition that is achieved (see Section [3.2.2t . yielding somewhat lower r], espe- 
cially for the innermost radial Lagrangian bin. This simulation also shows 
the same "energy divergence" toward < in the outer parts at early times 
as in Figure |6] 

evolution toward a state in which high-mass stars have lower veloc- 
ity dispersions than low-mass stars (i.e., 77 > 0). Interestingly, we 
have found a few instances in our simulations where this does not 
hold true. Instead, there is what one might call "energy divergence", 
where part of system evolves toward a state where high-mass stars 
have higher velocity dispersions than low-mass stars (i.e., 77 < 0). 

We find such energy divergence only in some of our simula- 
tions, and only for the outer ~ 20 — 30% of the particles during 
the first few initial half-mass relaxation times. The r]{t, r) profiles 
shown in both Figures |6] and |7] provide examples of this. The en- 
ergy divergence appears transient, and it develops over the same 
timescale over which maximum energy equipartition develops in 
the core. Because of this, we interpret the energy divergence as a 
result of dynamical interactions in the core. These interactions eject 
to the outer regions of the system particles that have partially ther- 
malized at a higher velocity dispersion, and that tend to be heavier 
than those native in the outer regions because of central mass seg- 
regation. 

Support for this scenario comes from the fact that the two sim- 
ulations that clearly show this feature are the IMBH run (Figure|7} 
and a run that forms a BH-BH binary and starts from a high core 
density (Wo = 7; Figure|6j. In both cases, it is expected that main- 
sequence stars in the core are scattered out nearly independently 
of their mass jGill et al.ll2008l) . After several half-mass relaxation 
times, this transient disappears once partial equipartition has time 
to develop in the outer regions as well. 

We also find < in some of our simulations at very late 
times, when the system has lost more than 80% of its mass. How- 



© 2013 RAS, MNRAS OOO.fTlfTTI 



No Energy Equipartition in Globular Clusters 9 



ever, this may very well be because very few main-sequence stars 
are remaining, so that the measurement of 77 suffers from signif- 
icant numerical noise. We therefore do not attach much physical 
significance to those instances of energy divergence. 



4 APPLICATION TO OMEGA CEN 

With the advent of HST proper motion studies of the internal kine- 
matics in globular clusters (see Section[T), it is now becoming pos- 
sible to actually measure the amount of equipartition in globular 
clusters. With the models presented in this paper, this opens up 
the possibility of detailed data-model comparisons. As a first ap- 
plication, we consider here the case of the galactic globular cluster 
Om ega Cen (NGC 5139). 

[Anderson & van der Marell (1201 0() determined a proper mo- 
tion catalog for the core of Omega Cen from ACSAVFC data with a 
4-year time baseline. The relation between a and m for stars along 
the main sequence with m = 0.5-0.8Mq is shown in their Fig- 
ure 25c. They showed that a power-law a = with rj ~ 0.2 
describes the data reasonably well. This implies that Omega Cen 
is between the extremes of no equipartition and complete equipar- 
tition. This is qualitatively similar to what we find generically in 
our A^-bo d y sim ulations. It is also consistent with earlier work by 
lAndersonI l2002h . which showed that there is mass segregation in 
Omega Cen, but not as much as predicted by multi-mass Michie- 
King models that assume complete equipartition. 

Andrea B ellini (priv. comm. 2013) r e fined the proper mo- 
tion catalog of I Anderson & van der MmcII (l2010t) in the context 
of a n ongoing HST Ar chival study of two dozen globular clus- 
ters ( lBellinietal1l2013h . He added more recent WFC3/UVIS data 
in many filters. This increases the time-baseline of the Omega 
Cen data to 8 years, and reduces all proper motion errors cor- 
respondi ngly. For the same sample of st ars as in the "central 
field" of [Anderson & van der Mare il (I2010I) , including stars down 
to m « O.3M0 but using improved proper motions, he finds that 
77 — 0.16 ± 0.05. This is the value that we will adopt here for our 
data-model comparison. The error bar includes the contributions 
of both random and systematic errors. Random errors are obtained 
from a least-squares fit to the a{m) data, with propagation of the 
random uncertainties for individual mass bins. Systematic errors re- 
flect, e.g., variations in best-fit rj values depending on exactly which 
WFC3/UVIS filters are included in the construction of the catalog. 

For a quantitative data-model co mparison, we need to estimate 
the dy namical age of Omega Cen. [McLaughlin & van der Marell 
l l2005h showed that the surface brightness is well fitted by a King 
model with Wo = 6.2l°;^, c = 1.311^3. ^ Ul."2+%'f, 
and Lv.tot = 10.0^-^°^° °^ . Ivan der Marel & Anderson! ilOldt) 
derived from a dynamical model for the spatially resolved kine- 
matics a distance D = 4.7 ± 0.06 kpc and mass-to-light ratio 
M/Lv = 2.64 ± 0.03 in solar units. From these values we in- 
fer using equation ([TJ a relaxation time trh = lO'''^'' yr at the 
half-light radius rh = 6.8 pc. For an assumed age of 13 Gyr, this 
implies t — 1.43 trh- 

For comparison to our simulations, we need to know the age 
in units of the initial half-mass relaxation time trh(0), and not the 
current half-mass relaxation time trh - Because of a combination of 
tidal evaporation of particles and contraction/expansion of the sys- 
tem, there is not a well defined relation between current and initial 
relaxation time. In the case of a system similar to Omega Cen, it is 
possible that trhif) is either larger or smaller than trh{Q), primar- 
ily depending on the retention fraction of stellar mass black holes 
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Figure 8. Data-model comparison of energy equipartition in Omega Cen. 
The time evolution of rj (the power-law slope in the relation between veloc- 
ity dispersion and mass, a oc m~'') is shown as in Figures|3]|6]and|2] Data 
points show the results from the A'^-body simulation snapshots for the inner- 
most 10% of the main sequence stars in projection. All A^-body simulations 
are included, and are generally shown in red; the simulation that includes 
an IMBH is shown in blue, and the simulation that is similar to that but 
without the IMBH is shown in magenta. The observed r} for Omega Cen 
is plotted with error bars in green at the current dynamical age t/trh{0)- 
Overall, there is excellent consistency between simulations and observa- 
tions, supporting the view that globular clusters are not generally in energy 
equipartition. Data-model comparisons of rj have the potential to help con- 
strain globular cluster dynamics, stnjcture, evolution, initial conditions, and 
the possible presence of IMBHs. 

(and on the presence/absence of a central IMBH), and on the ini- 
tial concentration of the system. If either stellar mass black holes 
or an IMBH are present in the system, then likely there has been 
expansion of the half-light radius fueled by energy generation in 
the core of the system, which implies trh(t) > trh(Q)- On the 
other hand, if Omega Cen lacks massive dark remnants, it is pos- 
sible that its half-light radius has undergone contraction since the 
cluster formation. Conservatively, we assign an uncertainty of 35% 
to the dynamic age of the system: t = (1.43 ± 0.5) trhiO)- This 
uncertainty includes both contributions propagated from the cur- 
rent model parameters that describe Omega Cen, and from the sys- 
tematic unknowns on the composition and expansion/contraction 
history of the cluster. 

To compare the observed amount of energy equipartition in 
the central field of Anderson & van der Marel (2010) to the simu- 
lations, we use the simulated Lagrangian bin that contains the inner 
10% of the particles in projection. When the simulation is scaled to 
the size of Omega Cen, these particles have a similar median pro- 
jected radius as the observed fieldjfl The comparison between the 

^ This comparison does not account for all the details of the observations 
and of the cluster For example, the field observed for proper motions is 
not a circle. Omega Cen is elongated and rotating, unlike the simulated 
clusters. And finally, the mass function in our simulations was not tailored 
to fit specific observations of Omega Cen. 
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measured r) for Omega Cen and the predictions of all of the TV-body 
simulations is shown in Figure [80 

Given the dynamical age of Omega Cen, the simulations pre- 
dict r] to be in the range 0.07-0.15. The measured r) = 0.16 ± 0.05 
falls within this range. So overall, there is excellent consistency be- 
tween simulations and observations. The simulation with an IMBH 
predicts a slightly lower rj than observed, while the simulation 
with the same initial conditions but no IMBH fits better How- 
ever, a broader suite of simulations with IMBHs would be needed 
to place any quantitative constraints on the possible presence of 
an I MBH in Omega Cen (a topic that continues to be d ebated, 
e.g.. lNovoIa et alj ( l201(]hfvan der Marel & AndersonldlOIOh ). With 
more sophisticated future analyses, smaller observational uncer- 
tainties, or larger samples of clusters, it may become possible to 
use observed rj values to discriminate between different clusters 
models. Detailed data-model comparisons of equipartition there- 
fore hold the promise of opening up a new observational window on 
globular cluster dynamics, structure, evolution, initial conditions, 
and the possible presence of IMBHs. 



5 DISCUSSION AND CONCLUSIONS 

It is widely believed, and commonly taught, that a globular cluster 
evolves, given a sufficiently long time, toward a state in which its 
stars are in energy equipartition (at least near the center, where the 
relaxation times are shortest). If the mean kinetic energy (imn^) 
becomes independent of mass m due to two-body relaxation (i.e., 
collisions), the velocity dispersion a = (i;^)"'^ scales as ct oc 
m""'^. Some popular multi-mass dynamical models for globular 
clusters, the so-called Michie-King models, have this scaling built 
in by assumption (Gunn & Griffin 1979). 

We have shown here that this paradigm is incorrect. The lumi- 
nous stars in a globular that has evolved for a long time converge 
toward a oc m"'', with 7700 ~ 0.08 ± 0.02 independent of position 
in the cluster. In this state, the velocity dispersion is nearly indepen- 
dent of mass. The inner regions can reach values up to r^max ~ 0.2 
after several initial half-mass relaxation times. Either way, the lu- 
minous stars in globular clusters with realistic IMFs always remain 
far from complete equipartition. 

The physical mechanism that explains why some syst ems 
cannot attain energy equipartition, namely the ISpitzeil ( [T969I) m- 
stability for a two-compone nt system, has been know for a long 
time. Also, Ivishniaj jl978l) showed that for a continuous mass 
function of the Salpeter form, energy equipartition is not attain- 
able. However, the implications of these results have gotten lit- 
tle attention in the literature. This is probably because previous 
investigations focused solely on th e stability of s ystems with re- 
spect to e nergy equipartition (e.g., ISDitzeilll969l [vishniadfigT^ : 
iKondrat'e v & Ozemov 1982) or the implied mass segregation. 
The latter is a consequence of energy equiparti t ion but depends 
on other things as w ell (e.g. iKhalisi et alj|2007l ; [oill et alj|2008l : 
IPasQuato et al]|2009l) . 

No previous theoretical study appears to have addressed the 
dynamical evolution of realistic star cluster models with a specific 
focus on mass-dependent kinematics, to determine exactly how 
close or far one expects systems to be from energy equipartition. 
We therefore made this the topic of the present paper This analysis 

^ Most of our simulations have initial concentration Wq = 5 or 7, which 
is close to the current value Wo = 6.2 j|g' j for Omega Cen. 



is particularly timely, because with the advent of HST proper mo- 
tion studies, it has now become possibl e to actually measure the re- 
lation between cr and m in real clusters dAnderson & van der Marell 
I2OIOI) . Measurements of this relation fo r main-sequence st ars in 
some two-dozen clusters are forthcoming ("Bellini et al .l20I3l) . This 
opens up a whole new discovery space. Detailed data-model com- 
parisons of mass-dependent kinematics have the potential to shed 
new light on globular cluster dynamics, structure, evolution, initial 
conditions, and the possible presence of IMBHs at their centers. 

We analyzed the stellar dy namics in the dire ct N-body simu- 
lations previously carried out bv lTrentietaI.I(l20i0>) . We quantified 
the relation between a and m, as function of time and position in 
the cluster, for realistic IMFs and initial conditions. We find that 
this relation is well-fit by a power-law of the form a oc , both 
for single main-sequence stars and for compact remnants. Com- 
pact remnants tend to have higher r/ than main-sequence stars (but 
still Tj < 0.5), due to their steeper (evolved) mass function. In 
the present paper we have focused mostly on the main-sequence 
stars, which are actually observable. Our main conclusions from 
this analysis are as follows: 

(1) The value of generally increases linearly in the first few ini- 
tial half-mass relaxation times (trh{0))- The increase is faster for 
Lagrangian bins closer to the center, where relaxation times are 
shorter. The central bin (containing the inner 10% of the stars as 
seen in projection) reaches a maximum ?7max ~ 0.15 ± 0.03. 
The increase to this value is mostly completed by f « 3tr/i(0). 
At large times, all radial bins convergence on an asymptotic value 
r/oo ~ 0.08 ± 0.02. 

(2) No simulated system ever reaches a state close to complete 
equipartition with 77 — 0.5. Even our most favorable conditions 
for equipartition to develop in the core, that is a steep IMF and low 
retention fraction of remnants, still yield only r^max ^ 0.19. Also, 
restricting the analysis to particles at the center of the system in a 
three dimensional, rather than projected, sense does not change our 
conclusions (although this does yield a slight increase Ari ~ 0.05). 

(3) Partial energy equipartition develops in an overall strikingly 
similar way across all our numerical experiments, despite the va- 
riety of initial conditions employed (compare, e.g.. Figures [3] [6] 
andlTJ. The maximum and asymptotic ri values do not differ much 
between runs (e.g., see Table[T]and Figs. |4] and |5j. Some trends are 
present depending on IMF and the content of compact remnants 
and binaries, and those trends can be understood based on simple 
physical arguments (see Section |3.2.2t . 

(4) The simulation with a central IMBH has the least amount of 
equipartition (as measured by r/max) among the sample of initial 
conditions considered. This result is consistent with the suppres- 
sion of mass segregation that we have observed in sim ulations with 
a cen tral IMBH iTrenti et al. 2007; Gill et al. 2008; P asguato et al] 
|2009). Further investigations with a larger sample of runs (espe- 
cially with an IMBH) are required to fully characterize the gen- 
erality of this result. Either way, this does suggest a new method 
for constraining the possible presence of an IMBH in a globular 
cluster, namely through the slope r] of the a-m relationy This is 

^ Searching for an equipartition signature from an IMBH in the cluster 
core has the potential to provide a diagnostic at earlier times in the life of a 
system than does mass segregation, since mass segregation is a consequence 
of equipartition. In fact, t > 5 <r?i(0) is required to discriminate s ystems 
with or without an IMBH from mass segregation l lGiU et alj|200g) . while 
potentially a signature in rj can be seen already at t > 2 trh(O). 



© 2013 RAS, MNRAS 000,[T1QT] 



No Energy Equipartition in Globular Clusters 1 1 



completely independent fro m methods based on the radial vel ocity 
dispersion profile cr(r) (e.g. Ivan der Marel & Andersonll201ol) . but 
can be pursued with similar (proper motion) datasets. 

(5) Any existing results derived from dynamical modeling that as- 
sumed complete energy equipartition by construction may be af- 
fected by unknown biases that will need to be carefully evaluated. 
For example, past studies have often relied on multi-mass M ichie- 
King models (e.g. Jde Marchi et alJl200d ; lBeccari et alj|201(]|) . and 
the underlying equipartition assumption of such models does not 
appear to be correct. It would be worthwhile for future studies to 
quantify any biases that might be introduced, or to modify the un- 
derlying model assumption to cr oc m^^ with rj 7^ 0.5. 

(6) Comparison of our simulations to a measurement of r) from 
HST proper motions in the core of Omega Cen yields good agree- 
ment. This supports the view that globular clusters are not generally 
in energy equipartition. 
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